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Preface 


We have decided to compile a preliminary documentation of the new 
GLAS Fourth Order General Circulation Model. The present documenta- 
tion has not been subjected to a careful editing process; we hope 
that its possible usefulness will compensate for some of its defects. 
The model dynamics (COMPO, C0MP1 and COMP2) is still undergoing 
minor improvements, especially in the time differencing scheme which 
we hope will improve its efficiency. The "physics" routine (C0MP3 ) 
has not been documented because it is being thoroughly revised. The 
present version of C0MP3 , similar to the one used in the GLAS/GISS 
models (see the documentation by Tsang and Karn), with modifications 
introduced by Y. Sud (.1979) is included in the code. Criticisms and 
suggestions for improvement will be reatly appreciated, since a final 
documentation will be prepared in 1980. 

We are very grateful to all the people that have helped us generously. 
In particular. Dr. N. Rushfield had a major impact in the process of 
making the model operational. W. Connelly, D. Edelmann, D. Han, 

S. Breining, P. Anolick and M. Almeida were very helpful, in the devel- 
opment of the model. The documentation was expertly and cheerfully 
typed by S. Mathis; D. Edelmann and D. Rosen have also cooperated in 
its compilation. We want to express our special gratitude toward 
Dr. Y. Sud, who offered us generously both his advice and his time 
in the development of the "physics" routine, and most especially to 
Or, M, Halem without whose many useful suggestions, constant encour- 
agement, and lpng patience, this work would not have been finished. 


November 1979 


Eugenia Kalnay-Rivas 
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I. Introduction 

The band fourth order model is a GCM which uses guadraiically conservative, 
fourth order horizontal space differences on an unstaggered grid and seoond 
order vertical space differences with a Matsuno (forward-backward) 
or a smooth leapfrog time scheme to solve the primitive equations 
of motion. 

This program numerically solves these equations one latitude hand 
atatiro'j which greatly reduces the amount of computer core storage 
needed to run the program. It also uses the same variable names, 
order of computations, I/O, post-processing as the standard second 
order GCM. Appropriate modifications have been made for the fourth 
order differences and leapfrog scheme. (See the 1978 Goddard Model- 
ing and Simulation Research Review for an overview of the fourth 
order band model.) 

The main feature of this model is that fourth order approximations 

gq 

are used for all the horizontal derivatives. The derivative z — is 

u X 

approximated by 

4 ,q (x+ Ax) -q (x-Ax) v _ 1 ,q (x+2Ax) -q (x-2Ax) , 

3 1 ZAx ' “ 3 ' 4 Ax > 

3 

and the derivative (qT) by 

dX 

4 r (T(x)+T (x+Ax) ) (q (x) +q (x+Ax)> - (T (x) +T (x-Ax))(q (x) +q (x-Ax)), 

j l 4^x 1 

1 r (T (x) +T (x+2Ax))(q (x) +q (x+2Ax)) - (T (x) +T (x-2 Ax))(q (x) 4-q (x-2 Ax)) , 
3 1 8 Ax J 
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The second approximation is derived by averaging the flux qT to 
yield a conservative form of the dynamic equations. Note that if 
T is equal to 1 the second equation reduces to the first. 


The primary variables are the horizontal components of the wind 
velocity, \V=(u,v), the temperature, T, the specific humidity, q, and 
the shifted surface pressure, tt, (7T=P s -P t0 p, P top = *° m b) • 


The secondary variables are the geopotential ,<p , the vertical wind 
velocity,©, and the pressure, p. 

The following ,pages give the differential equations of motion for 
the GCM model with the initial and boundary conditions. This is 
followed by the equations with the corresponding fourth order approx- 
imations which use the same notation as the current second order 
model. A complete description of the primitive equations with 

the o coordinate system is found in the Arakawa UCLA notes (1976) . 
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II . Primitive Equations of Motion 


1.&2. Horizontal momentum equations 


+ V. uw) + (ttcw) 


* -frV<p - Tra^I V7T~ (f+u^2l) IkxirV + 7r(F 
p a 


3. Continuity equation 

(3.1) + V • ( tiW) + ^-(ttcj) = 0 , or 


(3.2) 


T- -( v * 


(Ti\V)dcT = -V • C Ti\Vda 

'0 t) 


4 . Equation of state 
RT 


a = 


5. First law of thermodynamics 


+ V. ( rt VT) + |^(7rST) = m + 


P 

From 0 =T/p k , p=p T +cnr, u>=cnr+7ra , Vrr , k 

9 ttcfT k 9 TTO0 . _aair _ . c 

Ta~ = p TcT + 7T— (J~ Replacing in 5, 


9 ttT , , _k 9 7TCF0 _ 7T akT ,9 tt , ltr . ttQ 

3t + V(li\VT) +p -jj- = — — + \VVlr) + g* 

h* 


6. Humidity equation 
^ + V- (TiWq) = 0 


7. Hydrostatic equation 


9jj>_ 


9 k 
P 


-C 0 
P 


(from 


LE = -p = - 

9 cf> p 


=R/C 

P 


iH | <3 


Of the variables 1T,u,v ^ .T,q,<j>,a,o , we update the 5 primary variables 
•jr,u,v,T, and q using equations 1,2, 3. 2, 5 and 6. From equations 
3,1, 4, and 7 we can obtain 4>,ct, and o which are our secondary 
variables. Note that p=cnr+p t0 p. 

Sea level pressure (used only in the smoothing routine SMSHAP) 

Hydrostatic eg. * - ‘P - “ - -§f & - - f $ 

<P " S 

SLP *a p(o==l) exp{-~) 

RT 
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DERIVATION OF THE EQUATIONS AT THE POLES 


Consider the continuity equation 


It ♦ v * (7TV) + Itt ' 0 


dT 


coupled with a conservation equation jj^=S which can be expanded into 


9 nT + V. (WT) + 2*511 =71 s 


at 


3a 


If we integrate this equation over a polar cap of radius A<J> 


2/7 A <t> 


2 7T 


f ii? 
j 3 1 ** 



a^cos^difidA = - | ( V.(7iWT)a cos<j>d(j>dA 

0 J 

2tf tt^2-A4> 2tt jt/2 - 

SttcT - 2 COS(J)d( j )dA J L* -2 

/ 7T?2 



3a 


f r 

a^cos^dcJidA-} /7TSa & cos<|>d<f>dA 


9 ttT 

and we assume the value of to be approximately constant over the 


polar cap 


2 7T tt/2-A<J> 2tt 7t/2-A$ 

f j IF a^costftdipdX ~ (§^| f \ 

0 n/2 0 tt/2 


a cos<J>d<f>dA 


= 2?ra 2 (l-cosA<}>) (|^-) NP 

The first term in the rhs is, using Gauss" theorem 

2tt tt/2— A 4> 2tt 

-f ^ V. (ti^/T) a 2 cosc}>d<})dA = -J'iwT a sinA^dA 
0 ir/2 0 
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This can be approximated as 
IM 

~ra ,, i v i T i a sinA * 

The third term on the rhs is 


2 7T tr/2-i 


a 2 cos<J>d<f>dA - 27ra 2 (l~cosA*) (^~) NP 


and similarly with the source term. 


From these equations we obtain 


/3ttT v _ / i>m cotA4>/2 V* , m x 

‘rr’poie = - ( - 1) - alM ‘"iVi’ir 

1=1 T"^9 


. / 97raT v q 

( Ta~'Pole + S 



since 


2na sin A<(/2,r a 2 U-cosA*> = 2s i-nA »/2cos A^/2 = cot^| 
IM / IMa 2sin^ ^ alM 


This formulation is used in the continuity, momentum temperature and 
moisture equations. Note that the first term changes sign in the 
South Pole (m=l) . In the momentum equations we make use of the 
transformation 

IK = -sinAu^ - (-l) m cosAv^ 

VK = (-l) m cosAu^ - sinAv^ 

where m=l for the South Pole and m=2 for the North Pole. IK, \K are 
the "cartesian" velocities at longitude A^, and u^, the corres- 
ponding spherical velocities. 


The pressure gradient terms are computed making use of Green's 
Theorem : jj|£ dxdy =J’ Q d y; j'j|| dxdy = - Jpdx 
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For example 


2tt tt/2-A<J) 

f I 

I 5^ a 2 cos<J>d<J>cU * -$> <f>aA<J> (-cosAdX) 
0 tt/2 

4> m y — A$ 

IM 


+ lT & 


In the U-momentum equation we have then the following pressure terms 

IM 


>3ttUv _ A<J> 

'3t ; P a <i-cos4))XM f— . -'i 

i— x 

and similarly for ttV . 


£ + TiJcosXi 


P 'P i 


In the model we have approximated A<j>— 2sin2^ 


Then A(|> 

a (l-cos4>) IM 


a sin IM 


Based on this formulation we construct the fourth order scheme at 

4 

the Poles by taking j of the differences evaluated at A<£ from the 
Poles (as expanded here) , minus j of the differences at 2A<J> from 
the Poles. 

This formulation is not conservative at the Poles. However we have 
found that this has had no noticeable effect in the conservation of 
mass or energy in the model. In our shallow water experiments we 
studied a set of equations that were quadratically conservative, 
but inconsistent at the Pole, and another scheme analagous to the 
GFDL scheme, which is both quadratically conservative and consistent 
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at the Poles, but suffers from a serious truncation errcr near the 
Poles in the pressure gradient term. The scheme that we chose gave 
better results than ft’ae other two (Kalnay-Rivas , 1976). 
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Computation of the horizontal pressure gradient as suggested by 
N. h. Phillips 

(1) Let 6 * IT 4 O', F «= 280°K/1000 k i.e. constant 

(2) <$> * 4> 4 <£ ' , k * R/Cp * .286 

* ~Cp0 a 0 " <t> Q - CpFp k with q * CpFl000 k 

(3) T * T (p) 4 T' , T (p) * Fp k 

Thus our new dependent variables are 
<f>' = <t> + CpF (p k -1000 k ) 

T' * T - 0p k 

o’RT * 

In this way tt ( V4> 4 Vir) , the pressure gradient in the momentum 

r 

equations gets transformed into 

!c 

Tf ( V J 4 VJ>' 4 aReE- Vn 4 

J V P P 

_ 

= TT<V^' 4 SfiSlvff) 4 tt ( Vji" 4 aRQ^— Vu) 

But thfj (i?V,'ond parenthesis is zero: 

k , _ k 

Vi" 4 aR0 v tt— -Cp0Vp K 4 Re£-Vp = 0 

In regions of steep orography, the second parenthesis in (4) is much 
larger than the first. When the horizontal pressure gradient terms 
are computed in their original forru, the near cancellation of the 
two terms introduces large truncation errors. The procedure sugges- 
ted by Phillips greatly reduces this truncation error. We have chosen 
a simpler definition of <f> than the one suggested by Phillips. 
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III. 


Finite Difference Variables and Grid 


The notation used in fourth band model is the same as the standard 
GLAS(GISS) second order GCM except that we use a non-staggered hori 
zontal grid. A complete description of the variables can be found 
in the TSANG-KARN documentation of the GISS 9 level model. 


Let u ijk- v ijk’ T tjk- q ijk' *ij be the 

finite difference approximations to the primary variables u,v,T,q, 
and Tt at the mesh point (iAA,;)^* (jc^Ad) . Also the scaled variables 
ttu,ttv, ...are approximated by tiu^^, irv^^' ••• • 

The finite difference equations also use the following geometric arrays? 

DXP(j) = nij = a cos<f> AA, DYU(j) » n^ « aA<|K , 

Vie use a factor of 12 in DXYP to make our scaled fourth qpder differences simpler . 

Dxvp(j) = u-i-’j-nj , n A . = Dxyp(j)-,I.. , u. jk = n^u.^ 

®ijk “ n ij ,4 ijk ' ADLDP - 12 - a * 4XA * ' v Ijk " m j"ij v i 3 k 
P. k - DXYP(j)f. + ADi-DP-SIN*. U iJtj 

Horizontal Grid 


The fourth order baud model uses an unstaggered grid in the hori- 
zontal direction. 
j=JMl (N . Pole] 


A* £ 

j=l (S . Pole)^ 




" 








**" i 















L u "" 


_j 




. 




. 

=1 

i=2 

A A 

i=] 



JM intervals between 
south pole and north pole 


IM intervals in the longitudinal direction and periodic, 
i^IM+l is identical to 1=1. 
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A * = ^ 

(i-1) 5°-18Q° = ( j~l) A4i - ~ = (j-l)4°-90° 

JM1 - JM+1 

Vertical Grid 

The vertical grid is staggered; the values of all. the variables 
u,v,T,q,ir, ... except o are computed at the center of each layer. 
The values of 6 are computed at the edges cf the layers. 


a\ = llo:; * ~ 

IM 

A i ?= (i-1) AA-tt = 
JSP = 1 JNP = 



In the case of uniform vertical resolution and NLAY=9 vertical layers, 
a k = SIG(K) = (0 , 1/9 , . . . 1) 

k' = h>IGE (K) = is = 
and A o = 1/9 . 


j k and its scaled version S i; . k are defined at the eight interior 

edges, i.e. for o k with k=l to 8 since 6(0)=6(1)=0 from the 

boundary conditions. The pressure P k t sr '™ k ,+P TOp is defined at the 
same level as o k ,_^. 

Note for level k we need S. and S. . , , to form the second order 

13X X3/C-1 

vertical differences: 


! ijk ^ijk^ijk+l* - ^jk- 1 (V i j k +V ijk-l ) 


ORIGINAL PAGE «S 
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I I 1 - 2 Periodic Filtering of Short Waves 


An integral part of the numerical scheme is the periodic applica- 
tion (every ISMTH steps, generally 2 hours) of a 16th order Shapiro 
filter. This has the effect of removing waves shorter than 4 Ax, 
which are not resolved in the model, while waves longer than 4Ax, 
which are accurately computed by the difference scheme, are not 
affected by the filter (Kalnay-Rivas and Hoitsma, 1979). 

The filter is applied to an array q^ in the following way: 

Let d, (q . ) = q . , , -q . , d (q . ) ~ q,-q. , 

+ 3 4 J + i *3 ~ *3 J l“l 


Then a Shapiro filter of order 2N is given by 


N 


N 


q . = q . - (-1 ) (d, d) q. 

J J 

The response of the filter applied to a wave of the form 

q.=Q exp(i-Axj) is 
3 L 


,, . 2N ’ir Ax , 

(1 - Sin -jr~) q^ 


The 2-dimensional filter is applied as a product of 1-dimensional 
filters (first in longitude, then in latitude). In latitude we 
filter the fields on great circles formed by meridians of longitude 
A and A+ir, where 0< A <tt. We are presently filtering only potential 
temperature and sea level pressure. These fields were chosen because 
they are not very affected by orography. Winds are not currently 
filtered, because the adjustment between mass and velocity fields 
does not allow the development of short waves in the winds alone. 
However, in the tropics, where the adjustment of winds to the mass 
field it minimal, the winds are somewhat noisy, and we may opt to 
filter them too. 
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IV . Boundary Conditions 


Periodicity in the Zonal (East-We st) direction 

^ IM+m j — n mj n,=s ' j=2 , . . . , JM 

®IM+mjk~®mjk in^l #2,3/.../ j=2 , . . . , JM , k— 1 , . . . , NLAY 

for Q=U/ v,T,q,0,<5r, ttu, ttv, . . . 

Boundary conditions at the north and south po les. Define the array 
INDEX as follows: 



Then we can define 

7r iJM+2 = 71 INDEX (i) JM for t '" ie posits needed "beyond" the North Pole# 

and = 7T INDEX(i)2 for the P° ints "beyonc " the South Pole 

and similarly for T,q,<}>. 

For the horizontal velocity V=(u,v) we have 
V iJM+2k = ~ V iJMk 



V. 


Finite Difference Equations 


THE ZONAL (U) MOMENTUM EQUATION 


9 TTU 

3 t 


1 r 3 (7ru»u) + 3 (ttvcosi}) * u) , 3 (nau) 

acos <£ 1 3"X ’ §lj> J ” 3 a 

tt r 3 <f> , aRT 3 tt , . , j, . u tan <t> k 

5555* fjT + — JT 1 + <f + S~ ±)l,v + ffP x 


PU 1 i + l /2 


arrti COMP1: HA: PU1. .... 

^ u iik r 1-1/2 i - ** 

k}— = ( 4.*I (UjTT + U* , .. ) (U. .. +U. , .. ) - (0* +U?7,">t ) (U. . u +u.^, .. 

3 t ^ 13k i-l^k 13k l-ljk 13k i+l3k 13k i+l3k 


)] 


PU2 


i +1 


'• 5 1 <U i jk +U i- 2 jk> (U i jk +U i- 2 jk> ‘ <U i jk + ^i+ 2 jk> (U ijk +U i + 2 j k> 


PV.1 


+4f 1 ,v Ijk +v lj-lk ) (u i j k +u ij-lk ) ' (v Sjk +v l j+lk ) <°ijk +u ij+lk )1 


PV2 


• s *'< v *i:k +v * j . 2k ) (u ijk + 0 ij- 2 k , - <v Ijk +v !j+ 2 k> <" ijk + » 1 J+ 2 k> 


VA: 


+ '5* [ ®ijk-l <U ijk +U ijk-l )_ ^i j k <u i j k +U ijk+l ) ) / Ao k"l 

COMP 2 : PC: a RT? . 

+ (tt . . *n ,* { 8 * [<j> ! n ..-4> ! ,, .. +— — 2 Jji ^i-lj _ 7 T i+l j ^ + 
13 3 r i-l 3 k r i+l 3 k p, J J 


d. RT! 'I 

[ *i + 2jk-n-2jk + -Tf- (T i + 2r’'i-2j )1 ) 


C: 


COMP 3 


+ F. . ,*TT. . *V. ) + dlF ) 

13 k 13 13 k nc X 


Note: COMP1, COMP 2 , COMP3 are the names of the three subroutines 

where the different terms are computed. HA: Horizontal ad- 

vection terms. VA: Vertical advection terms . PG: Pressure 

gradient terms. C: Coriolis term. 

Also note that PV2 is set equal to zero for j=0 and j=JM, i.e., 
there is no transport of mass across the poles. 
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THE MERIDIONAL (V) MOMENTUM EQUATION 


3itV 

at 


1 r a ( 7i u • v) , a (ttvcos4)V , 3tt6v 

~ 'V r* a i “ g g 


acos4> 1 ax 


3 <j> 


~ * — It 1 - (f+ — TTU + 7fF 

a 3 4> p 3f a y 


COMPl: HA; 


anv. .. 

llii 

at 


{ 4 -*‘< U Ij k +U I- 1 jk) (Vijk+Vl _ 1Jk )-(OI Jk +Oi +ljk) < Vljk +V 1+lJk ) 


-• 5[(U *jk +U *i-2jk ) * (V ijk +V i-2jk>-' U ijk +U * + 2jk> (V ijk +V i + 2 j k ) 
+4 • * 1 (V * jk +V *i j-lk» (V ijk +V ij-lk) - <Vt jk +V^ +lk ) < V i j k +V i j + lK ) 


5* [ (V* jk+ VJ j+2k ) (V ijk+Vi j. 2k ) - (VJ jk+ Vj j+2k ) (V . jk+ V ij+2k ) 


VA: 

+ . S* [S ijk .i (V A jk+ v ijk . 1 ) -S . jk (V . jk+ V . Jk+1 ) ] /Aa k } 


COMP 2 ; P: 


0, RT! .. 


+ U i .*m.M8.*[4.! j . lk -*! j+lk + — p^- k - 


U;,IU , ,, 

+ [c ^i j+2k”^i j"2k ^ pT 3- (7T i j+2 =1T i j-2 ) ] } 


C: COMP 3 ; 

■ F ijk**ij* U ijk )+ <*F y > 
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THE THERMODYNAMIC ENERGY EQUATION 


3 irT 
3 t 


1 r 3 ttUT 3 ttVcos^ »T i k 3 ttctT/p^ 

acos<f> 1 3 \ 3 4 > a a 


ttoKT ,3ir 


U 


3tt , V 3tt 


11 + ^ on , v uir. ttQ 

v 3t acos (J) 3X a 34>' + C 


OOMP1: HAs 


anT, 


^ = { 4? [ (U^ +UJ . .. ) (T. . , +T. . ) - (UJ . +U* i (T +T \ 

3t *3* i-ljk' i 3 k i-ljk' ^ u ijk^ u i+ljk M1 ijk +T i+ljk ) 

-.5*[(U* +U* ... ) (T, ,,+T, - (u* +U* ) IT +T \ 

i]k i-2jk i 3 k i-2;jk' 1 ijk i+2jk' 1 i jk +T i+2 jk' 


+ 4* t(VJ jk+ Vj j . lk )(T.. k+ T. j . lk ) - <VJ jk +VJ. +lk ) <T ijk +T. j+lk ) 


• 5 [ (V ijk +V i ;) _2k l <T ij k+T ij+2k ) " <V ijk +V ij+2k ) <T i jk +T ij+2k > 1 
VA: 

p ijk p ijk-l p ijk p i jk+1 


ijk 
COMP 2 : 


q 

^ijk J 


+ V *jk {8 -* [ ’'i j+ l^i J -l] + ’ r i j - 2 -’'i j+ 2 }I > 


COMP 3: 
nQ. ■ 


+ (■ 


±jJS 

Cp 
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THE MOISTURE EQUATION 


g-ig = 

3 t 


9 nq. 

3 t 


_i_ ( ^g + iiv co af a , . *i2a + n!E . c) 

COMP1 : HA 

(4.*[(UJ. k +U*_i jk ) (qi jk + qi.i jk )-(uj jk+ uj +ljk )(q. . k+q . +ljk ,] 

-.5*((UJ jk+ UJ +2jk ) (q iJk +q i _ 2j!t )-(OJ jjk+ UJ +2; . k ) (q ijk+ q i+2jk ).l 

+4 -* I<V !jk+ V I j -lk ) << >ijk + « i j. 1 k>-<VI jk WI j+1It ) (q ijk +q ij + lk n 
-• 5 * I (v! j k + vi j . 2k > (< 3ijk + '3i j -2 k > - <v ij k + vj j+2k ) (q ijk + q ij+2k > ) 

+ - 5 * ‘®i jk-1 (q i jk^i jk-1 1 -Si jk ‘li jk+^i jk + i> 1 / Ao k } 

COMP 3: 

+n(E ijk- c ijk> 


Note: The current transport scheme for the moisture field is being 

modified. 
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THE PRESSURE-TENDENCY EQUATION 


3 TT 

at 


r * ft r airU x iiveosd), 
A=1 acos<f> 3l 3<f> J 


COMP1 : 


C0NV iftj =A^{8.MUj. ljr UJ +lj ,) + Uj +2jr Uj. 2j , 
+ 8 - MV ij-ir V ij+U )+V ij+ 2 r V ij-2A } 


COMP1 : 
L 
l 


a n. . 

■f? ■ l til C0NV Uj 


THE VERTICAL VELOCITY EQUATION 


air 1 r dnV t dnVcostp, 3 y „ ix 

ft = ' TX + ■ diT 1 * So (,T0) 


Thus 

A0 k = CONV ikj - AS ijk 


giving: 


COMP! : 


• • g jj 

S ijk = {S ijk-l + C0NV ikj " Aa k 8T } 
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VI. 


The Forecast Equations at the Poles 


The fourth order band model uses a spherical cap at the poles, and 
the finite difference approximations to the equations of motion must 
be derived for this spherical region. Stereographic projection is 
used to give us a well-defined velocity vector at the poles. 

From trigonometry, (Fig. Via) if the vector (U,V) in X-Y coordinates 
is represented by (U',V’) in X'-Y' coordinates where the prime co- 
ordinate axes are rotated by an angle X then, 

U) U * U’cosA - V’sinA 

(2) V = U’sinX - V'cosA 

We want to interpret “polar stereographic" 
velocities (U Np ,V Np ) also as "spherical" 

velocities (U,V)^, as we approach the 

Fig. via North Pole along a meridian of longitude X. 

From Fig. VIb, which shows a unit vector in both coordinate systems, we 
see that, at the North Pole, spherical coordinates are rotated by an 
angle (ff+A) with respect to the polar stereographic coordinates. Therefore, 
7 

U Np as Ucos(A+j) - Vsin(A+j) 

V Np = Us in ( A+ j) + Vcos(A+~) 

or 

(3) 

Similarly, for the South Pole: 

(4) 

also, 
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Fig. VIb 








The initial values for U Npk , V fJpk , U gpk , V gpk 
are obtained from equations 3 and 4 by averaging in the zonal 
direction on the line of latitude j*45 for the north pole and j=2 
for the south pole. In the forecast stereographic velocities axe advected 
and transformed back into spherical velocities after each time step 
by equations 5 and 6. 

In the following equation we will denote the polar velocities by 

VP^ m , VP^ m and the temperature and other variables in a similar 

way, TP. . ttP. , ... where m=l for the south pole and m=2 for the 
iv 9 in k t in 

north pole. The following constants are used, 

COEFl=(-l) m C0EF2=-C0EF1 RADIM=3*a.iM 

C0N1--4C0T ( . 5A4>) /RADIM 

CON 2 *■ -COT (A <J> ) /RADIM 

CON3=4DT/(RADIM*SIN( . 5A4> ) 

C0N4=DT/ (RADIM* SIN (A<p) ) 

JPOL (K ,M) - f 2 JM ^ 

U JM-lJ 

r=jpoL (1 ,m) irst interior value of j ( 2 for the S. Pole, JM for tne N. Pole) 
s=JPOL (2 ,m) =second interior value of j(3 for the S. Pole, JM-1 for the N. Pole). 
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ZONAL (U) MOMENTUM EQUATION (POLES) 


9 7TUP, COMPl : IM 

- { C0EPl*IC0Nl* i E i -7,. r V. rk (U irk SinA i+ C0EPl*V il . k C0SX i : 


IM 


+ CO N 2 *^-\ s V isk ( U lsk Si -' A i + COEFl * V lak COSA i )] 
1 - 1 

+ . 5 [fiP k . 1 (UP k «P k . 1 ) -SP k (UP k +UP k+1 ) ] /4c k ) 


COMP 2 j 0 pipp * 

+ ^m 1 ^ *if x ( n rk +JS pP^ *ir> C0SX i 


+CON 4 * 1 ^! < A i °k RTP k ff.JCOSAJ 
2 — — — is 1 


DT 


1=1 


isk pP. 


COMP 3 : 

+ f m 7 TP m * VP ^ rj + ( 7TP m ‘ FP V m 

mm k , :n m k , m 


MERIDIONAL (V) MOMENTUM EQUATION (POLES) 


COMPl : 


IM 


3 ttVP 

k,m ={COEFl* [CONI* £ ^ ir ^ irk (COEF2*U irk *COSA i +V irk SinA,.) 

a t i = o , 

IM 

+ CON2*J i - nis V isk (COEP2*U isk *COSA. + V isk SinX.)J 


+ . 5 [SP k _ 1 (VP k+ VP k . x ) -SP k (VP k+ VP k+ 1 ) ! /A 0 k ) 

C 0 MP 2 : jj, a RTP 

. r CON 3 A “ /A , , k KIF k , 

+ I Krn-*^ ^i r k 71 ' -) SmA 


m 1 DT 


i=l 


IM 


+ con 4*^ (<i>; 0l , +c7kRTPk 


DT ^ 1 ' r isk pP^ is 
C0MP3 : 

“ f ttP *UP, m + ( ttP * SP, ) 
mm k ,m m k/in 


P p k ir 
•tt . )SinA,] 
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THE THERMODYNAMICS ENERGY EQUATION (POLES) 


COMP1 : 


3 nTP 

TF 


k,m 


IM 


IM 


" tCOEFl*(CONl* * » lr V lrk T lrk *001l2» I 


TP TP 


pp 


k,m pP k-l,m 


TP. m TP. .. m 
-SP,_ _< — )J/Ao. ) 


r “TT 

pP k,*n pP k+l ,m 


COMP 2 : 

TlP 0 , KTP, „r3TTP_ 
m 

P P ! 


TlP 0 KTP. ( 3 ttP 

+ ( — SLJS — + CON 5* E I (COEFl*UP,. .COSA.+VP 


k,m 

COMP 3 : 

ttP 

+ <c p -o p k> 


IM 

E 

i=l 


k, m 

( 8 ’ T ir-” is ) ) ) > 


i k, 


THE MOISTURE BALANCE EQUATION (POLES) 


a COMPl : TM 

3 Tiq.P. IM 

= (COEF1* (CONI* E 


IM 


3 1 


i=l 


Tr ir V irk c 


irk 


■fCON2*E 

i=l 


7T is V isk 


^isk^ 


+ . 5 [SP k-1 (q p k , m +( 3 p k-l ,m ) ~® P k ,m (qP k ,m +qP k+l ,m ) ] /Aa k J 


COMP 3 : 

+ (7T V EP k ) 


SINA . )* 
m i 
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THE PRESSURE TENDENCY EQUATION (POLES) 


XM 


IM 

CONVPL v =COEFl* (CONI* E tt. V. .+CON2* E tt 4 „ V, .)Ao. 

XT lrK J-S 18 K K 


3 7TP f 

IT 


m 


NLAY 
E CONVPL 
k=l 


k ,va 


SP 


THE VERTICAL VELOCITY EQUATION (POLES) 

3-rrP 


k , m 


SP k-l,m + C0NVPL k * Ac, k 3 1 


rn 
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VII. Diagnostic Equations (<P,a,p) 


Once the updated values of itU, nV, nT , irq are found we unscale: 


..n+1 ,. T n+l . n+1 
U ijk * ffU ijk /ir ij 


for all i,j,k. 


Similarly for V,T,q. We also filter the fields near the poles to 
prevent linear instability (see subroutine AVRX) . 6 is obtained 

from S by unsealing also. 

n-M n +l 

We determine from and 


n+1 _ _ n+1 

p ijk - Vij + p 


TOP 


P TOP =C ° nStant * 


(ju is the surface geopotential (a function only of latitude and 
b longitude) 

k 

For the Phillips geopptential we define (p) at the center of the layer in the 
following way: 


k+1 k+1 

(?) k , P ijk+r " P ijkV „ 

Ajk (k ‘ fl) (p ijk'+l^ p ijk^ 


(P 


1 1 n, k + 1 

k_ 1 3 p ^ 


k+1 9 p 

k=l , 2 , . . . , NLAY 


where 


" <SKE(k)* I i lj+ P Top ) k+1 


(SIGE (k) =a£) , i.e. 


k+i . 


is obtained by exponentiation and (p) bv differences, 


The following equations represent the geopotential calculations 
used in the old fourth order model. 
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Lot 


C 

ijk 


w ii°k RAq k 


p i jk 

for k«l, . . . , NLAY 


C 

2 ^ 


k 

v’ / p i jk-fl 
k'l k “ 
p ijk 


l) + o k .. x (l 4 ^i 

/ ' Pi-ik 


ijk / J 


with 


k - k k 


P ijO p ijl and p i jNLAY-f X = p 


ijNLAY 


An optimized version of C. is* 

13k 


c 

ijk 


n ij °k RA0 k _ ^ 5 C E 


ijk 


P 


ijk 




^ 0 k' (l, ijk + i-Pijk )+ “k-i<p ijk+ i-p- ;jk . 1 )| 


Rather than compete * and then subtract ¥ wc do everything at once 

U) ^ijNLAY = *' B - CPTHMPSKAPA - p* ) + T’c . . . T . . „ 

1JLLAY ^ 13 d ij£ 


<2> *ij» 


where 
CPTM = C 


0 


^ij Ji-tl + (p ij e-^l'“ p ±jk 


PSKAPA « 1000 k = p k 


) ( + ^itl) _ 20 , 

P ijk P i jk+1 
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VIII . 


The Time Differencing Scheme 


The fourth order band model has the option of using the Matsuno 
time scheme or the smooth leapfrog scheme (see MWR-Vol 100 (487- 
490) R. Asselin) . 

Let Q n represent a typical variable that is to be updated to time 
n+1, and let D(Q n ) represent the nonlinear space differences. The 
Matsuno (Euler-backward) scheme is as follows: 

Q = Q n +AtD(Q n ) 

Q n+1 « Q n +AtD(Q) 

The standard leapfrog scheme is given by 


( 1 ) 


Q n+l_ Q n-! 
~TSt 


D(Q n ) 


For the smooth leapfrog scheme we replace Q n ”"^ by Q n- ^ 

(2) Q n+1 =Q n “ 1 + 2 AtD (Q n ) 
with 

(3) Q n = (l-v)Q n + .5v(Q n ” 1 +Q n+1 ) 

Equation (3) represents a simple time filter except Q is used 
1 

instead of Q r in order to save core storage. The above equa- 
tions (2) and (3) represent the order in which the smooth leapfrog 
scheme is evaluated. For n=l we define Q^=Q° then we update in 
equation (2) followed by the filtering in equation (3) which is 
needed for the next time step. 
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The smoothing step introduces dissipation with respect to time, 
controlled explicitly by the parameter v, as compared to the im- 
plicit dissipation in the Matsuno scheme. The amplification 
factor can be found in the paper by Asselin. 

A further modification must be made to the smooth leapfrog scheme 
when source terms are included. Essentially the idea is that we 
must include the source term effect over two steps rather than one 
If we do not do this, then the source effects (C0MP3) are included 
only in every other step which will introduce large discretization 
errors. (For details see the attached report. Appendix B.) 

If the source terms are called every NC0MP3 steps, then for 
step n=NC0MP3 

Oi - Q n ~ 2 + 2 AtD (Q n_1 ) 
q£“ 1= (1-v) Q 11 ” 1 + .5v(Q n " 2 + Q n ) 

Then compute the source terms £> n ^nd include in both steps n and n 

Q n = Q* + NC0MP3 * AtS n 
Q n “ 1 = Q*” 1 + NC0MP3 • AtS n 

The actual code is complicated by the fact that we actually use 
scaled and unsealed variables, but the generalization is straight- 
forward. 
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IX. Documentation of the Code (Preliminary) 


The band fourth model uses special equivalences and nonstandard dimen- 
sions in order to have the variables P,U,V,T,SH stored in contiguous 
locations for each line of latitude. Thus, we desire to have the 
variables stored as follows. 


IM NLAY*IM IM NLAY*IM IM 

» 1 1 — * 1 »- 

1 IM+1 (NLAY+1) *IM+1 
PHIS . u . TS , v SHS 

<— — * 1 *- 

The scaled variables PT, UT, VT, TT, 


IM NLAV*JM IM NLAY*IM IM 

> f » » * 

1 IM (NLAY+1) *IMfl 

PU . UT , GT . VT , GW , 

— 1 1 4 — — + 


NLAY*IM IM NLAY*IM IM 
* , 1 h 




SH 


•PHIS . U 

•— jif- ' " 


SHT are to be stored as follows: 
NLAY*IM IM NLAX*m IM 

, I » 1 


TT 


PT 


SHT 


PU . UT ... 

h -jzi ' 


The above storage designation is accomplished by dimensioning 
PHIS (4* (NLAY+1) *IM,JNP) instead of PHIS(IM,JNP) (for a fine grid 
we have PHIS (2880 ,46) , for ultrafine we have PHIS (4800 , 72) ) . Then 
we equivalence NLAY*IM locations of the first line of latitude of U 
with PHIS (IM+1, 1) to PHIS ( (NLAY+1) *IM,1) (the first IM locations of 
PHIS are used for the quantity PHIS itself) . 

Similarly, we equivalence TS(1,1) to TS(IM,1) with PHIS ( (NLAY+1 ) * 
IM+1,1) to PHIS ( (NLAY+2) *IM, 1) and so forth. (See the enclosed com- 
puter codes for the exact values in the fine and ultrafine versions.) 
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In order to have the successive lines of latitude arranged properly 
in storage we dimension our variables to achieve this purpose. Thus, 
we have U ( IM, 4* (NLAY+1) , 1) instead of U (IM, NLAY, JM) , similarly for 
V,T,SH,UT,VT,TT, and SHT. For P we have P (4* (NLAY+1) *IM, 1) instead 
of P (IM, JNP) and similarly for TS, SHS,GT,GW, and PT. Note that the 
variables U,UT,... are only computed at IM*NLAY*JNP points, the 
special dimensions are needed to properly align the variables in 
storage. It is important to recall that v/e are using two proper- 
ties of the FORTRAN Compiler: 

First, by equivalencing two elements of two different arrays we im- 
plicitly equivalence the other elements of the arrays. Second, that 
last array dimension can be left as 1 as long as it is dimensioned 
properly in the calling routine or it is equivalenced to a properly 
dimensioned array. Note we could have dimensioned U as 
U (IM,4* (NLAY+1) , JNP) we would use the same amount of storage. It is 
crucial to have the IM* (NLAY+1) dimension because it would cause the 
computations on U for the successive lines of latitude to be shifted 
IM* (NLAY+1) locations where the next line of latitude of U are stor- 
ed. For clarity and simplicity in programming we use the standard 
equivalence of U,V,T,SH with Q(I,L,N,J) with U equivalent to Q(I,L,1,J) 
and so forth. Since we want P,U,V,T, and SH in contiguous storage 
locations and the equivalence of Q with U,V,T, and SH, then we must 
dimension Q as Q (IM, NLAY+1 , 4 , 1) instead of Q ( IM , NLAY ,4,1) . We fill 
in the extra locations by including PHIS,TS,SH, and P. 
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COMP0 Description 


This subroutine contains the time schemes and controls the callina 
sequences of the routines COMPl, COMP 2 , COMP 3 , and the polar filter- 
ing routine AVRX. 

The logic of this subroutine involves three considerations . First, 
it permits one to choose either the Matsuno or the smooth leapfrog 
time scheme. Second, the latter scheme involves the use of storage 
arrays PSM and QSM. Third, the unsealing, smoothing of the updated 
variables QT and the call of subroutine COMP 3 occur at the value 
JS2=J-2 when J is the value of the current line of latitude which 
is being computed. The routine COMPl and COMP 2 require the values 
of U , V , . . . at JS2 in order to compute the updated values at J because 
fourth order differences are used in the meridional (J) direction. 
Only after the COMPl and C0MP2 are called for value J can we unscale, 
smooth, and finish processing the variables at JS2. 

The code for the poles is identical in format with the code for the 
other J values except the variables are scaled by the pressure PPOL 
only . 

The main program contains two calls to COMPJ^ which cannot be treated 
independently because of the calling sequence: CO^P0(Q,QT) f COMP0 

(QT,Q) . If the first call is a leapfrog step (LF) , the second one 
must be also LF. If the first call is Matsuno predictor (MP) , the 
second one must be Matsuno corrector (MC) . This is represented 
symbolically by LF-^LF, or MP->MC. The second call can be followed 
by either MP or LF. Each of these combinations requires different 
transfers. 
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Description of the Time Step Sequence Parameters 


NSTEP: Counts the time steps. Starts and restarts both begin 

with NSTEP=0 . 

NSEQ: The number of steps (combined matsuno and leapfrog) in 

each (repeated) sequence of time steps. 

MLF(I): MLF ( I ) = 0 or 1 according to whether the Ith step in the 

sequence is Leapfrog or Matsuno, respectively. First step 
is always Matsuno (MLF(1)=1). 

ISMTH , \ Smoothing routine (SMSHAP) is called MOD (NSTEP-NSM1 , 

NSMl: | ISMTH); if XSMTH=0 , there is no smoothing. 

NCOMP 3 ,\ Physics routine (COMP 3) is called MOD (NSTEP-NCM1, 

NCM1 : j NCOMP 3 ) . 


Sample Runs 
Matsuno only: 

NSEQ*1, MLF (1 ) =1 , MATSUN=1 , DT=750., NSM1=0, NCM1=0 
BCIN03=4, ISMTH=8 


Leapfrog only: 

NSEQ=1 , MLF(1)=1, MATSUN=0, DT=600., NSM1=0, NCMl=0 
NCOMP 3= 5 , ISMTH=10 

1 Matsuno, 4 Leap-Frog: 

NSEQ=5, MLF= (1,0, 0,0,0) , MATSUN=not needed, NSM1=0, 
NCM1=0, DT=600. , NCOMP3=5, ISMTH=10 
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C0MP1 Description 


The COMP1 subroutine contains the horizontal and vertical advection 
differences. The DO loops over I are arranged to make use of the 
periodicity of the variables in the zonal (I) direction. For example, 
suppose we are to compute D(I)=Q(I*fl)-Q(I) for 1=1,..., IM. Then the 
corresponding code is 

I=IM 

DO 10 IP1=1 , IM 
D(I)=Q(IP1)-Q(I) 

I = IP1 

10 CONTINUE 


Where we used the periodicity Q ( IM+1) =Q (1 ) . 


In the meridional (J) direction we compute our difference approxima- 
tions in stages in order to make maximum use of each line of latitude 
of a typical variable when it is in core. The fourth order differ- 

g 0 

ence approximation to ( i£ j A<f>) =D (i , £ , j ) 

(1) 4 /^iftj+1 ®i£j-l v 1 ,^i£j+2 ^i£j-2^ 

- j( L_) 


Thus we see that Q, g . will be needed in the difference approximations 

J 

9 O 

to for j-2, j-1, j+1, and j+2. The corresponding code is 

DO 20 j p 1 IM 

QFLUX1= 4* (Q(I,L,JP1)+Q(I,L, J) ) 

QFLUX2=- . 5* (Q(I,L, JP2)+Q(I,L, J) ) 

D ( I , L , JP2 ) = D ( I , L , CJP2 ) -QFLUX2 
D ( I , L , JPl ) = D ( I , L , JP1 ) -QFLUXl 
D (I , L, J) = D { I , L , J) +QFLUX1+QFLUX2 
20 CONTINUE 

For simplicity the array D is inititalized to zero and Q is scaled 
so that (1) contains no divisions. 


At the poles (j=l or J=JNP=46) we use the values given in section IV 
on boundary conditions. We have special code for these cases denoted 
J=2 or J=JM corrections. 
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COMP 2 Description 


This routine contains the Coriolis force term, the geopotential 
calculation (which should be made into a separate subroutine) , and 
the pressure gradient and energy term calculations. 

The geopotential PHI is dimensioned PHI (72,9, 5) since we only need 
at most five storage locations for any computation. We use modular 
arithmetic (M0D5) to compute the indices JMOD, JP1M0D, JP2M0D which 
correspond to the standard index values of J,J+1, and J+2. Thus 
PHI (I, L, 6) is stored in PKI (I,L,1), and we avoid shifting array values 
by using the JMOD index as a pointer. 

For the south pole calculation we need geopotential values at J=2 
and 3, thus, for J=1 we compute PHI for J=l,2, and 3. For successive 
values of J we need only compute PHI at JP2 which is needed in the 
pressure gradient calcualtion at J. Therefore, the calculation of 
PHI and the associated array PK are coded for calculation at JP2. 
Except for the first J value we are only computing the geopotential 
at one latitude value for each pressure gradient calculation. 
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X: 


Flow Charts 


CGMPO (Q,QT) 


JS1=1 

JS2=1 


( (MAIN LOOP) ) 

DO 10 J=1,JM 

( (COMPUTE ALL J-EARAMETERS) ) 


IF (JP2>JM) 
IF (J>1) 


GO TO 25 
GO TO 18 


JP2=JP2MQD=2 

( (SAVE QTPOL (M) IN QSMPOL (M) ,M=1 , 2) ) 
18 ((SAVE QJ(JP2) IN QSM ( JP2M0D) ) ) 


IF (JP2>2) 

( (INCREMENT JP2, JP2MOD) ) 

25 CALL CQMPl (Q,QT,J) 

CALL COMP2 (Q,QT,J) 

( (ELIMINATE NEGATIVE HUMIDITIES) 

IF(P'l4(400,1100.)) 

IF (J<3) 

IF (J=3) 


GO TO 25 
GO TO 18 


STOP 

GO TO 200 
GO TO 70 


29 ((UNSCALE QT(JS2))) 

CALL AVRX (QT, JS2) 

( (LEAP FROG, MATSUNO AS RELATED TO THE 
SEQUENCE OF 2 CALLS TO COMPO: 

1st CALL ^ LF $ LF or MP =* MO, 

2nd CALL ^ LF 4 LF or LF 4 MP or MC ^ MP) 

IF (LF^LFj GO TO 45 

IF (LF^MP.OR.MC^MP) GO TO 58 

( (STATUS: MP+MC OR MC+LF) ) 

( (P(I,JS2)<-PSM(I, JS2MOD) ,Q<-QSM) ) 

IF (MC+LF) GO TO 63 

(STATUS: MP+MC) ) 


X-l 


(contin. ) 



COMFJ (Q,QT) 


58 

45 

63 


64 


( (P(I,JS2HPT(I, JS2) , Q»-QT*DX¥P*PT) ) GO TO 64 

( (P(I,JS2)<- *P(I,JS2)+a* (PSMfPT) , 
Q^3*DKVP*P*Ofa* (QSM+QT*DXYP*PT) ) ) 

CONTINUE 

( (SOURCE TERM CORRECTION FOR LEAP PROG) ) 


IF (NOT C0MP3 CALL) 

( (0<O-OT*DXYP*PT) ) 
IF (NOT C0MP3 CALL) 

IF(MATSUNO PREDICTOR STEP) 


CALL CQMP3 (QT,JS2) 

( (COMPLETE SOURCE TERM CORRECTION, 
(M>K2T*DXYP*PT) ) 


GO TO 67 

GO TO 67 
GO TO 67 


67 IF ( J<JM) 


GO TO 200 


( (INCREMENT JS2,JS2MOD) ) 

IF ( JS2<JM) 


70 

((POLES)) 

200 

JS2=JS1 


JS1=J 

10 

CONTINUE 


((END CF J LOOP)) 


FEIURN 


GO TO 29 


COMF1 (Q,QT,J) 

( (COMPUTE JP2-JS2M0D) ) 

IF(J>2) GO TO 2150 

JS1=JS2-JS1MQD=JS2MQD=1 


(contin . ) 
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C3QMP1 (Q,QT,J) 


2150 IF (J=OM) 

((COMPUTE PV W PV 2 =V* + VV 0 )) 

% 3 3 +* 

IF(J=1) 

2158 ( (COMPUTE PU^, PU 2 = U* + U| +2 > > 

((COMPUTE HORI Z . ADVEC . - IN IONG. DIREC. 

IF (J = JM) 


2225 

2290 

2237 

2405 


( (COMPUTE HORIZ. ADV.-IN LAT. DIREC.)) 

IF <J>1) 

( ( J=2 CORREC. TO HORIZ, ADVEC & CONV. ) ! 
((CONV. CALC. FOR CONT. EQ.)) 

IF (J<JM) 

( (J=JM CORREC. TO HORIZ. ^VEC. & CONV', 
• 

((CONV OONV * DSIG(L) ) ) 

IF (L<NLAY) 

IF (J=.l) 

((COMPUTE SIGDOT,PT) ) 

( (COMPUTE VERTICAL ADVEC. ) ) 


IF (J<JM) 


2600 


((POLES, M=1 or 2)) 
((HORIZ. ADVEC.)) 
((CONT. EwQ. ) ) 

( (SIGDOT PL, PTPOL) ) 
((VERT. ADVEC.) ) 


END OF COMP1 


Alternative code 

from 

to 


■ • — 

— 

IF(J.EQ.l) 

( (COMPUTE PU] 

i, pu.n 


((COMPUTE HA' 

(i) ')) 

2222 

( ( J=2 CORREC 

)) 


2235 IF ( J < JM) 

( (J=JM CORREC) ) 


2225 ((HORIZ ADVEC (J) )) 

( (CONV. CALC. FORCONT.EQ. ) ) 

2405 ( (CONV CONV*DSIG(L) ) ) 

# 
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GO TO 2158 
GO TO 2225 

) 

GO TO 2237 
GO TO 2290 
GO TO 2405 

GO TO 2150 
GO TO 2600 

RETURN 

RETURN 

GO TO 2222 
GO TO 2235 

GO TO 2225 
GO TO 2405 



COMP 2 (Q,QT,J) 


3001 

3005 

3007 

3055 

3065 

3030 

3032 

3065 

3135 

3111 

3031 


( (COMPUTE JP2=. . . ) ) 


IF ( J>1 ) 


( ( JP2=JP2MOD= s JPKP2=l) ) 


GO TO 3001 
: 3005 


( (CORIOLIS ) ) 


IF(J_>JM) 


GO TO 3032 


( (FIRST MAIN LOOP IN L) ) 


IF (L<NLAY) 


( (DO 3030 LX=1 ,NLAY ) ) 

( (COMPUTE ENERGY TERM IN 
THERMODYNAMIC EQN ) ) 

( (L=NLAY+1-LX ) ) 


GO TO 3055 


( (COMPUTE PK ( I , LL , JPKP2 ) FORLL= 1 , NLAY ) ) 
((COMPUTE PHI' (I,NLAY, JP2MOD) ) :3065 


((COMPUTE PHI' AT JP2MOD FOR L<NLAY) ) 


IF ( JP2=1 or JP2=JNP) GO TO 3030 

((COMPUTE W ( JP2MOD) ) ) 

( ( INCREMENT LX ) ) 

((IF J=1 GO TO 3005 AND COMPUTE PHI (2 ) , W (2 ) ) ) 


( ( 2nd MAIN LOOP IN L) ) 
((DO 3031 LX=1, NLAY)) 

IF(J=1) 

IF (J > 2) 

( ( J=2 s VT,TT OORREC.) ) 
(((P) FOR U HQ.)) 

IF (J=2orJM) 

( ( (P) FOR V EQ. ) ) 


IF (J < JM) 

(( CORREC. TO (P) J=JM) ) 
( (POLES) ) 

( (CONTINUE) ) 

((RETURN)) 

END 


GO TO 3111 ( POLES) 
GO TO 3085 

GO TO 3135 

GO TO 3031 
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A FOURTH-ORDER FORECASTING MODEL* 

(E. Kalnay-Rj vas, D. Hoitsma, and P. Anolick) 

The GISS fourth-order model { Kalnay-Rivas , et al . , 

1977), which is a fourth-order, energy-conserving GGM~on an 
unstaggered grid, had shown promising capabilities. It produced 
forecasts that showed an improvement over the second-order GISS 
forecasts with the same fine grid (4° x 5°) resolution, but that 
were somewhat inferior to the "uitrafi.ne" forecasts. However, the 
first version of the model required excessive amounts .of computer 
memory and time for execution.. 

The model has been reprogrammed into the "fourth-order 
band model." The new program solves the primitive equations one 
latitude band at a time. The arrays are stored in an interlaced 
way, with all arrays being updated at the same latitude stored 
contiguously, and similarly for all arrays used in the computation 
of the time derivatives. This design of the program makes effective 
use of the virtual memory capability of the Modeling and Simulation 
Facility's IBM 370/165 or Amdahl computers. The virtual memory 
facility permits the execution of programs whose cure size is larger 
than the one available, by placing the excess on disk and reading in 
those pages of information needed in the current calculations. The 
band fourth-order data structure and computations were constructed 
to optimize this virtual I/O process; a possible improvement may be 
to interlace also the arrays being updated with those used to compute 
the time derivatives. The use of the virtual memory facility avoids 
the explicit I/O used in the current Kern model, and yields a simpler 
program. (The band structure was suggested by G. Russell.) 

The band fourth-order model computations have been 
optimized so that each time step is computed in the Amdahl in half 
the time required by the old model. The array structure has also 
been designed to reduce the amount of overall storage and high- 
speed memory by a factor of two (see Table 1). 


Table 1. Comparison of 

Fourth-Order Model 

Computing Requirements 


Original 

B . tad 


4th-Order Model 

4th-Order Model 

Core (bytes) 

3 50 OK 

1500K 

CPU time per step 



(COMP1, COMP2) 

34 sec 

17 sec 


♦Reprint from Atmospheric and Oceanographic Research Review-1978, 
O r AS , NASA Technical Memorandum 80253, 
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The model has been programmed so that it can use both 
the Matsuno and the leapfrog time schemes, the latter with the 
Robert time smoother. If Q n and D(Q n ) denote the fields at time 
nAt and the corresponding time derivative computed from thd, space 
differences, the smoothed leapfrog scheme is 

Q n+1 _ gn-1 + 2 At D(Q n ) 

Q" - Q n + . SvtQ"” 1 + Q n+1 -2Q n ) 


with Q° « Q°. The use of a smoothing coefficient requires a 
slightly smaller time step. For example, with v = .1, the 
model is marginally stable with At = 288 sec. , compared to 
.Vi = 300 sec. for the Matsuno and leapfrog schemes. 

The smoothed leapfrog scheme is further modified to 
include the subgrid '’‘physics" terms, and scaling and spatial 
smoothing procedures. Since the "physics" is called every 
few time steps, unless the leapfrog scheme is restarted after 
every call to the "physics," only one of the two consecutive 
fields will be affected. The restarting procedure is time- 
consuming, so it has been replaced by the following algorithm 

Q n-1 . Q n-1 Q r. 

Q n . Q n + s n 

Q n-1 . Q n-1 + qU 


where S n corresponds to the "physics" terms. This procedure, 
which ensures that the physics is applied tc two consecutive 
time steps, has been tested with good results. 

The model can be extended into an "ultrafine" version 
in a straightforward way. 

Preliminary Results . The first numerical integrations performed 
with the new band fourth-order model show dramatic improvements 
over forecasts made with the GLAS model with the same resolution 
("fine grid" or 4° latitude by 5° longitude). The quality of the 
forecasts is now comparable with those produced with the "ultrafine 
(2,5° latitude by 3° longitude) version of the model. Results 
from a 3-day numerical integration are presented in Figure 1. 

During an extended 8-day integration of the new fourth-order 
model, the atmospheric systems remained remarkably smooth, 
exhibiting a realistic behavior both with respect to position 
and intensities. 
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Abstract: The new CilSS 4th order model of’ the global atmosphere is described, It is based on 4th order 
quudrutically conserving differences with the periodic application of a 16th order filter on the sea level pressure 
and potential temperature equations, a combination which is approximately enstrophy conserving, Several short 
range forecasts indicate a significant improvement over 2nd order forecasts with the same resolution (~ 400 km). 
However the 4th order forecasts arc somewhat inferior to 2nd order forecasts with double resolution, This is 
probably due to the presence of short waves in the range between 1000 km and 2000 km, which are computed 
more accurately by the 2nd order high resolution model, An operation count of the schemes indicates that 
with similar code optimization, the 4th order model will require approximately the same amount of computer 
time as the 2nd order model with the same resolution. 

It is estimated that the 4th order model with a grid size of 200 km provides enough accuracy to make 
horizontal truncation errors negligible over a period of a week for all synoptic scales (waves longer than 
1000 km), 


1. {introduction 

It is generally accepted that the use of 4th order finite differences is more efficient in reducing 
space truncation errors than the use of higher resolution on a 2nd order model (KRlilSS and OLIGI5R, 

1 072), Unear analyses and shallow water type of experiments give an upper limit of the ini|u ovement 
that can be expected from the reduction of errors in the horizontal differences. For example, Table I 
corresponds to a linear wave equation with phase speed e = 1 1 ms" 1 , typical of atmospheric motions, It 
provides a measure of the “computational predictability period” after which the errors introduced by 
horizontal truncation alone become very serious, The table suggests that for waves longer than 2000 km 
a 4th order 400 km grid model is preferable to a 2nd order 200 grid model. Waves shorter than 
2000 km are forecast more accurately by a 2nd order 200 km grid model. In order to insure that 
horizontal truncation errors are small in the 1-3 week period of atmospheric predictability for all 
synoptic scale waves it is necessary to use either 4th order differences with a grid resolution of the 
order of 200 km or 2nd order differences with a grid of the order of 100 km. Numerical experiments 
with simple nonlinear models (WILLIAMSON and BROWNING, 1973, KALNAY-RIVAS 1976a, from now 
on 1) also indicate dramatic reduction of errors by tire use of 4th order differences, 

All these studies indicate that a considerable improvement in forecasting skill is to lie expected from the 
use of 4th order differences if a substantial portion of the forecasting errors is due to horizontal trunca- 
tion errors In actual numerical forecasts there are several are several other important sources of errors: 
vertical truncation errors, errors in the initial conditions, and poor “physics” (i.c. parameterization of 
physical processes like radiation, dissipation, subgrid transports, cumulus convection, boundary layers, 
etc.) 

* (Paper presented during the DMG-AMS Meeting, Hamburg 1976: see Preface to Issue I 2/1977) 

3oitrd ge zur Physik dor Atm o sphiiro, 60. Band, 19?7,Seite 299-311 
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■ Table I , Flapped time T utter which u wave of wavelength L moving with a phase 
speed o w 1 1 ms" 1 lags by more than 100 km due to space truncation errors. 



Ax « 400 km 

Ax - 200 km 

2ml order 

t, * 2000 km, l 0,4 days 

L « 2000 km, T - t ,6 days 

differences 

l « 4000 km, 1 = (.6 days 

L • 4000 km, T ~ 6,4 days 

4th order 

1- 2000 km, T~ 1,5 days 

l~ 2000 km, T“ 21 days 

differences 

I.b 4000 km, T <- 21 days 

1.*= 4000 km, ,128 days 


At the Goddard Institute tor Space Studies (GISS), New York, we have developed a 4th order general 
circulation model (G(’M) with the expectation that it witl yield more accurate short range forecasts, and 
probably more realistic climate simulations, hut with the improvement limited by the other sources of 
errors. The model is described in Section 2, The results of several experimental short range forecasts are 
discussed in Section 3, A sample 36-hour forecast is presented and compared with forecasts by the 
standard 2nd order GISS model, which has the same resolution as the 4th order model, and by the 
“ultrafine” 2nd order GISS model, which has twice the horizontal resolution, Section 4 contains the 
conclusions and a discussion of future work. 


2. Description of the model 


The 4th order GISS global atmospheric model is a primitive equation model using longitude 
(X), latitude (0) and sigma (a) coordinates, The basic equations and the parameterizalions of physical 
processes are the same as those of the standard 2ml order GISS global model (Somerville et a!, 1674). 
The finite difference scheme is quite different and is described in the following subsections, 


2, 1 . Finite-difference scheme 

The horizontal grid is uniform (constant AX and A0) and non-staggered, There arc two ways 
to define such a grid, depending on whether variables are defined at the poles, as done by Wli 1 i AMSON 
and HKOWNiNG (1973), or half a grid size away from the poles (HOI I OWA Y et al, 1973), We chose the 
former method because in the* absence of smoothing near the poles it allows a time step twice as large 
as the latter. The singularity that shcrical coordinates have at the poles is explicitly avoided by the use 
of a polar cell where “stereographic” velocities are used. 

In the vertical direction we use a staggered grid with a , the vertical velocity, defined at the boundaries 
of each layer, and all other variables in the center of the layer, The vertical grid can be non uniform 
although experiments have been made sc far with constant An, (experiments reported in I indicated that 
no significant improvement in accuracy was obtained when a 4th order staggered conservative scheme 
was used in the vertical direction. Therefore the model has 2nd order vertical differences, The compu- 
tation of the geopotential <l> is performed as indicated by ARAKAWA (1 972), with a modification 
suggested by fill 1. 1,1 PS (1974). 


a) Forecast equation away from the poles 

Several systems of horizontal difference werl tested as reported in I, The scheme chosen for the 
b:»sis of these experiments consists of t lie simplest possible quadratically conservative 4th order differences. 


We define the finite difference horizontal divergence operator at the vertical level k as 


»k(g> = aJs0 [Jmh'uV) 3«2k(llu n g aA ) 


+ 1 5 ^ ( 1 1 v cos <;> 0 g 0 ) !. 8 2 \ ( H v cos 0 20 g 3 * ) j 

J k 


(2,1) 


and the finite difference gradient operator as 


t 7 ? ^ 

8 acos0 


3«2Ag 


2 A 


+ S|%g" ^2,g 20 ] 


( 2 . 2 ) 


Here II = ps pr is the difference between the pressure at the surface and the constant pressure at the 
top of the model, a is the radius of the earth, \v = uii + v g is the horizontal velocity vector in spherical 
coordinates, and we use the finite difference notation 

fi nAg = |g(X + n AX/2, 0, o, t) g< v~n AX/2, 0, a, t)|/(nAX), 

(2.3) 

g' ,X = |g(X + » AX/2, 0, o, t) + g(X n AX/2, 0, «, t)|/2, 


and similar formulas for the other independent variables. With this notation, the continuity equation away 
from the poles is 


an 

at 


D k (l) ll«„Ok 


(2.4a) 


or, if we integrate it in the vertical and make use of the boundary condition 
a “ 0 at 0 = 0, 1 , 

31 , K 

at “ V'MOAok (2.4b) 

k I 


where K is the number of vertical layers. Equation (2,4b) is the forecast equation for II. The momentum 
equation is 



, / iiif tan <6 \ 

D k («/) II a „ (0 V 0 \ + [ f + a ) IK X ( I hv) k 


II 



VII 


+ IF|c> 
k 


(2.5) 


The first law 6f thermodynamics is 

B - 0,.( T ) - {p* IW 0 \6 (T/p |} k 


+ 


IIokT 
P 


an 

at ~ 


+ w.VII 


4* 


(2.6) 
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Here T is absolute temperature, F is the frictional force, k ~ R/C,„ Q is the diabatic heating per unit 
mass, p is the pressure and p K is computed as indicated in Subsection 2,1c, The moisture equation for the 
water vapor mixing ratio q is 



l\lq> HM(Vq‘V ll(K C). 


(2,7) 


f and (' are the rates of evaporation and condensation, 

In the momentum equation (2,5) we follow a device suggested by mil I It'S (1074) to alleviate the 


difficulties of the computation of the pressure gradient in o-coordinates in tin* vicinity of orography . 

We define T' ~ T T(p), «l>' « <|> <|>(p), where T « 0 p\ fl = 280 K/(I000 mh) K , <I» - «1> 0 C„ Op'* and <h 0 = 
|( 1 000 mb) * pi- 1- \Ve have simplif led 1*1111 1 It'S’ expression for <1*0 since its precise value is not 
important. With this procedure, in the pressure gradient term there is an exact cancellation of the terms 

o K F 

V ‘I* + 11 VII, and this implies a significant reduction of truncation errors in regions with steep orography, 


b) 


Forecast equations at the Poles 


We have followed the method used by wit l IAMSON and imoWNlNC! (1073) and define a polar 
cap of radius a Ay> on which we use “stereogruphic" (or rather “cartesian") velocity components defined 
by the transformation 


Uii> » usinA t vcosX 
V , |* » ♦ u cos <\ vsin A 
with inverse 

u ~ U , |* sin * V,|) cos,\ 
v ~ • Hi,, cos A V) |> sin A 


(2,8) 


(2,‘>) 


where the top and bottom signs correspond to the north and south poles respectively. The positive 
x-axis of the cartesian coordinates coincides with the meridian of longitude A (A The difference in 
signs ai the south pole is due to the choice of a right handed system of coordinates with the verti- 
cal unit vector pointing outwards, formulas (2.8) and (23)) are used to define "spherical” or ‘‘cartesian’ 
velocity components wherever (hey are required in the finite difference equations, 

The finite difference horizontal divergence operator at the poles is defined by an average aver all 
longitudes Xj ; 


i i * 


!>>»’, k(‘?) ^ j *5 At (llvg),*.;, , 


iii/i) 


i a 


; A, (II vg)\j, , (’I 2A*)| 


( 2 . 10 ) 


where A| - (i 3) AA, AA s 2/r/l, and A„ - AAasin(n A0)/| 7 n a : ( I eos(n A0))|. The Unite difference 
gradient operator at i lie poles is also defined hy an average over all longitudes: 


1 I 2 


Vtpg : 


la A 0 


i a 


j , 8 a j, * (J M<) 6 8,\|, . (-J 2 m ,) j ‘ {cosA; : * sin Aj j] 


(2.11) 
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Then the continuity equation at the poles is 
K 


dlli j* » i 

»*^D t ,..k(DA<, k , 

k I 


( 2 . 12 ) 


the momentum equation is 

("”, V ) * * D t P,k ( V) - » . P*0 (o tp Vtp) + f±p to (II V) t (l , 
ll ll .[v t ,-‘l»'+“r V ‘I ,M ] +FiP.k 


*P“ U,tl 

i> i + V , |> j 

1. The 

(bl\T\ 



lit*) 

1 “ f D tl i )k 

tp, k 

til- 

0 K Ti|> 

■aiiti 

1 

Ptp 

af 


(2.13) 


(2.14) 


+ v±»* * v,p n 


ii*p 

+ Q±,, 

v n 


and the moisture equation is similar. 

c) Diagnostic equations 

In o-coordinates, the pressure is defined by p = oil + pp. Following a suggestion by PHILLIPS 
(1974) we define at the center of a layer p k = !/(« + 1) 8„p k + '/8 0 pk. instead ofp£ = (pk)* as in ARAKAWA 
(1972), This formula was derived under the assumption that the potential temperature varies in a o-layer 
much less than either the temperature of the pressure. PHILLIPS indicated that a more accurate relationship 
between temperature and geopotential can be expected from this formula. 

The geopotential 4> is determined following ARAKAWA (1972): If «l> s is the surface geopotential, then the 
geopotential at the center of the lowest layer is 


<1>k = 4‘a + VCkT k 

k “ I 

and at other levels 


/ T a \ 

4'k = 4> k + , + C„ Ao k + i fi 0 (pk + i ) ( I 

\ Pk + 1/2 / 

The coefficients C k are determined from 

llo k RAo k c„ — -r-a 

(«a o5oP ) k 


(2.15) 


(2.16) 


(2.17) 


The vertical velocity in o-coordinate is determined from (2.4a) and (2,4b), and the boundary condition 
b = 0 at o = 0, 1 : 


«0 = If { V D k (l) Ao k - D k (l) 

k = I 


(2,18) 
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d) Boundary conditions 

In ihe east-west direction we use periodicity: g(\ + 2 n,<j>) = g(X, <p) for all variables, When the 
value of a variable is needed “beyond” the poles, we define it by continuation along the same meridian: 

«(X, * ('i + &<P )}= V (i r + X ;g(x, + X* *(f “ Ao)) » 

where g represents all variables other than the two horizontal velocity components. In the vertical the 
top and bottom are material surfaces through which there is no ilux (6 = 0 at o = 0, I ) except for subgrid 
boundary layer fluxes of moment'^' and heat included in F and Q, 


2.2. Filtering near the poles and high order filtering 

The CFL computational stability condition requires a very small time step unless linearly 
unstable short waves are filtered out near the poles, For this purpose several alternative procedures 
have been tried but so far the method found to give most satisfactory results is the Fourier filtering of 
the prognostic variables u, v, T and the indirect smoothing of 1 1 through the filtering of the sea level 
pressure (SLP) field, 

The Courier filtering is performed polewards of 66° latitude. The amplitudes of the Courier components 
of zonal wavenumber n are multiplied by a transfer function which is I for n < N and decreases linearly 
to zero between n = N and n = N + 5. The number of retained modes is defined by N($) = integer part 
of (90 cos0). We have tried filtering the stereographic velocities U, V instead of u, v but uo improvement 
was obtained. 

Since both the surface geopotential and the surface pressure fields contain large amplitude short wave 
components due to the presence of orography, an artifical smoothing of these fields represents a 
distortion of the real geometry of the boundary, On the other hand, the SLP is an instrinsically smooth 
field, and its high wavenumber components more closely represent atmospheric waves. Therefore in our 
model the SLP is filtered near the poles and 11 is recovered by solving the transcendental equation used 
to relate them 

It was found that the 4th order model forecasts were less smooth than those of die 2nd order GISS 
model, especially in regions with steep orography. Based on the same considerations we have introduced 
a periodic application of a high (16th) order filter (SHAPIRO, 1970) on the SLP and potential tem- 
perature fields, which are not very affected by orography. 

The filter is of the formg = { 1 "(Fp 8 }{ 1 •~(F^) 8 }g, where F^(gij) “(gi + ij * 2g S j +gj-jj)/4 and has a response 
jexp(ikA)] = - sin 2 (kAX/2) exp(ik\). This filter eliminates waves shorter than 4AA and even after 
hundreds of applications has negligible damping effect on waves longer than 4AX, 

In the meridional direction there are three simple alternative forms for F^: 

F>i (gij) ~ (gij + i 2gjj + gjj . j )/4 

P? «Ml(Bj) = l(ftJ+ I ~ gij) COS0j + , /2 ~ (gjj “gij- i)cos^.., /2 1/(4 COS*j) 
and 

F 0lll (8ij) = I F 0l (gij cos^X/cos^j 

FRANCIS (1975) used P^ n , which is the only one that conserves the area weighted average of g, 

However, we have found that the three meridional filters produced virtually identical results, This could 


A-8 



be expected because the waves affected by the (liter are too short to be strongly influenced by the 
convergence of the meridians. Since Fj ( and Fj,n can be programmed ns efficiently as F* we find them 
preferable. In the model we use with g continued “beyond” the poles as indicated in Subsection 
lid. At the poles g t} . is defined as the average over all longitudes X| of the filtered values g*, *p. At the 
present we are applying the high order filter to the SLP and potential temperature fields once every two 
hours. 

2.3. Enstrophy constraint 

HAYLiss and ISAACSON (I97S) have developed a simple method to make any finite difference 
scheme conservative with respect to any quantity. The method has been tested in our model by forcing 
conservation of enstrophy on the dry adiabatic version of the model, although in such a model it is 
potential enstrophy that is conserved. The procedure is the following: Let the vectors U and V represent 
the values of the velocity components discretized over the grid and let the functional G(U, V) denote 
the consistent 4th order approximation of the mean square vorticity. At each time step a correction 
U\ V' added to the predicted values U, V such that G(U + l)\ V + V') = G(Uo, V 0 ) where U 0 , V a 
denote the velocities at time t = 0. Since this equation cannot be solved explicitly, it is linearized about 
the predicted values U, V: 


(iuW u, ^v)o. V-W*., V.)-0<O.V) 

(2.19) 

If we assume 


(U, V)-a( 3l j, 9 y )_ _ 

(2.20) 

then « can be determined from (2. 1 9): 


G(U<), Vo) ™G(U, V) 

1 i)G I 2 . | <)G | 2 

(2.21) 


oti * , tiu 

0 0 0 V 


It may be easily shown that the choice (2.20) minimizes II If 1 II + II V ' II .the norm of the correction vector. 
It has been found that this procedure improves the forecasting skill of the model. However, when it is 
used in combination with the periodic application of the high order filter the improvement is marginal. 
Therefore this option is not included in the standard version of our model. 


2.4. Summary of the properties of the model! 

The finite differences of the model have the following properties: 

a) Horizontal differences are performed on a nonstaggered grid and have 4th order truncation errors. 

b) Vertical differences are performed on a staggered grid and have 2;td order truncation errors. 

c) The differences are conservative in the sense that in the absence of diabatic and dissipation, mass and 
energy are conserved except for marginal terms at the poles. Non-conservative differences similar to 
those used by WILLIAMSON and BROWNING (1973), starting from real data, proved to be nonlinearly 
unstable after less than one day of integration. 

d) Unlike the “box method”, the horizontal differences remain 4th order near the poles (KALNAY-RIVAS, 
1976b) 



e) The model contains no horizontal diffusion expect tor the possible use of a dissipative time scheme 
and for the periodic application of the high order filter. It has been shown in I that the combination 
of a 4th order quadratically conservative scheme with the periodic application of a high order filter 
replaces successfully the use of an enstrophy-conserving scheme. This is because waves shorter than 
four times the grid size are the ones subject to aliasing and to large truncation errors, and they are 
removed by the filter before they attain finite amplitude. Waves longer than four times the grid size are 
accurately computed by the 4th order scheme and are not affected by the filter. 


3. Numerical forecasts 

3.1. Analytical initial conditions 

The model was tested using 3 dimensional extensions of both the nontrivial steady state 
solution and the Rossby-Haurwitz wave initial conditions used in I. The forecasts remained smooth 
during the several days of alteration. Errors in the steady state case were an order of magnitude larger 
than in I because of the use of single precision arithmetic on the IBM/360-95 computer. However, it 
was found necessary to use double precision for the longitudinally averaged pressure gradient at the 
poles in order to avoid excessive local error growth, 


3.1. Forecasts with real initial data 

As mentioned before, the parameterization of physical processes included in the terms F, 0, C 
and E of Equations (2.5) (2.7) has been adapted from the standard 2nd order GISS model (SOMERVI LEE 

et al, 1974). The only change that has been introduced is the reduction by a factor of 2 of the surface 

drag coefficient . This was found to be necessary in order to avoid an excessive damping of the pressure 
systems after one or two day forecasts. We believe that the fact that in the staggered grid the surface 
winds are obtained by a horizontal average reduces the relative effect of friction in the standard model. 

At the present time (November 1976), we have performed 5 experimental short range forecasts, Since 
they have been made with slightly different versions of the 4th order model, we don't yet have a reliable 
measure of the model’s forecasting skill, The initial data so far has been available in the staggered grid 
used by the standard model so that winds have been linearly interpolated to the nonstaggered grid. This 
is a source of error which may have adversely affected the 4th order forecasts. 

We have performed a single 4-day forecast which indicated that the model remains very stable and that 
synoptic systems show no tendency to become unrealistically weak, 

All the numerical 4th order forecasts show a significant improvement over the 2nd order forecasts 
performed with the same resolution. This improvement appears as a better estimation of the changes in 
position and intensity of several pressure systems. This in an encouraging result, especially in view of the 
study made by BAUMHEFNER and DOWNEY (1976) which indicates that the standard 2nd order GISS 
model forecasts compare favorably with those made by the NMC and NCAR models with similar 
resolution. 

On the other hand, and contrary to our expectations, the forecasts made with the 2nd order GISS model 
with double horizontal resolution were found to be either similar or superior to the 4th order forecasts. 

This is more true in the sea level pressure than in the 500 mb forecasts, We consider that there are 
several possible reasons for this result. One is the extra errors introduced in the 4th order model initial 
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• Figure I. Observed 36 hr, sea level pressure at 0000 GMT, 22 February 1 976 


data by the averaging of the winds. A second reason (and probably the most important one) is the 
existence of waves with significant amplitudes and wavelengths between 1000 and 2000 km. These 
short waves are better detected in the initial data and more accurately forecast by the 2nd order model 
with double resolution. The nonlinear interaction of these short waves with longer synoptic waves 
which according to linear theory are better forecasted by the 4th order model is clearly very important, 
A third important reason is tlw fact that errors introduced by the parameterization of subgrid processes 
become less important as the grid size is reduced, 

We present an example of a 36 hour forecast of the sea level pressure corresponding to February 22, 
1976, 1 2Z. Figure 1 is the verification sea level pressure map, Figure 2 is the forecast computed with 
the standard 2nd order GISS model, which has a resolution AX = 5°, A<j> = 4°, A o = 1/9. Figure 3 is the 
4th order forecast computed with the same resolution, and Figure 4 is the forecast obtained using the 
"ultrafine” 2nd order GISS model which has a resolution AX = 2.5°, A</> = 2°, Ao = 1/9. 

The most striking improvement showed by the 4th order forecast is in the position and intensity of the 
low in the eastern portion of North America, In fact, this low and the high over the Atlantic Ocean have 
been forecasted better by the 4th order model than by the 2nd order model with double resolution, 

On the other hand the developing cold high in central Canada and the intense low pressure center south 
of Greenland have been predicted better by the double resolution 2nd order model. As usual in this 


A-ll 




type of comparison, the three forecasts share several important deficiencies, For example the intensity 
of the high over southwest USA has been overpredicted by the three models, probably because of 
difficulties introduced by the orography. The low south of Alaska has been erroneously split by the 
three models. This splitting, which is slightly less pronounced in the 4th order forecast, may be due to 
both errors in the initial data and orographic and coastal problems, 


4. Summary and conclusions 

We have described the characteristics of the 4th order G1SS model of the global atmosphere. 

It is based on a quadratically conservative scheme with the periodic application of a 16th order filter on 
the sea level pressure and potential temperature fields, As shown in I this combination is approximately 
enstrophy-conserving. An operation count of the numerical schemes indicates that with similar code 
optimization the 4th order model will require approximately the same amount of computer time as the 
2nd order enstrophy-conserving GISS model with the same resolution. We also plan to introduce a 
simplified semi-implicit scheme, and to study the possibility of using the combination of a nonstaggered 
vertical grid with the Kreiss 4th order method described in 1. 
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The results of several short range forecasts indicate a significant improvement over the 2nd order forecast 
with the same resolution. This improvement is shown in the estimations of changes in position and 
intensity of several pressure rystems, We plan to study the impact that tire greater accuracy of 4th 
order differences has on the forecasting skill of variables of more practical importance, such as tempe- 
rature and precipitation. 

It lias been found that the 4th order forecasts are somewhat inferior to the forecasts made with a 2nd 
order model with double horizontal resolution. We consider that one of the most important reasons 
for this result are the presence of waves in the range between 1000 km and 2000 km which are com- 
puted more accurately by the high resolution 2nd order model than by the 4th order model. Another 
important reason is that errors introduced by the parameterization of subgrid processes become smalle r 
as the size is reduced. 

It is estimated that the 4th order model with a grid size of 200 km provides enough accuracy to make 
horizontal truncation errors negligible over a period of a week for all synoptic scales (waves longer tlu,n 
1000 km), 
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• Figure 4. 36 hr, using the “ultraline” 2nd order GISS model with double horizontal resolution 
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). INTRODUCTION 2. DISCUSSION OF THE DESIGN OF THE MODEL 


Tha design of a numerical model for ataoapharlc 
alaulatlon la not a straight-forward procedure. 
Both In the areas of oatheaatlcal and numerical 
analysis, and tn the parameterization of 
physical processed not explicitly resolved, 
the modeler faces several difficult choices 
between equally reasonable methods, and 
sometimes between similarly unsatisfactory 
■ethoda. 

In this paper we discuss the considerations 
leading to the numerical design of the CLAS 
Fourth-Order Global Atmospheric Model, This 
medal, which was briefly described in Kalnay- 
Rivas at al, . (1977], has been restructured, 
and several minor changes were Introduced, 
flit computation time and memory requirements 
for the 4 *h order model are now similar to 
those of the present second order GLAS 
model with the same 4° latitude , 5° longitude 
and 9 vertical-level resolution [Soimerville 
tt al . . 1974], However , the fourth-order 
model forecast skill is significantly better 
than that of the current GLAS model, and 
after 3 days it is arable or better than 
that obtained with the 2, 5° by 3° version 
of the GLAS model. 

A discussion of several of the basic charac- 
teristics of the model design Is contained 
In section 2. For each of them we present 
some of the possible alternatives, their 
advantages and disadvantages, and the reason 
for our choice. In section 3, we discuss 
the effect on numerical forecasts of changes 
In the accuracy, resolution and conservation 
properties of the model*.' 1 ' Section 4 
contains some final remarks. 


(1) As of this writing, several of the 

numerical experiments are not complete. 

Present affiliations Department of 
Mathematics, Wesleyan University of Connecticut 


Different numerical analysts and atmospheric 
modelers often take different approaches In 
the design of a numerical model for weather 
Prediction, Sometimes they even have d-fff*r*nt 
basic philosophies. For example Arakawa (1966, 
1972], has been a pitcneer In the development 
and use of numerical schemes that reproduce 
as closely as possible the conservation pro- 
perties of the continous equations of fluid 
dynamics that determine the motion of the 
atmosphere. On the other extreme, Krelss 
and Ollger (1972, 1973) have advocated the 
use of more accurate schemes even when they 
don’t formally satisfy any conservation 
properties. 

The current GLAS model Is based on Arakawa 'a 
(1972] sc iond order scheme using a staggered 
grid B. The scheme Is energy conserving and 
approximately enstrophy conserving for non- 
dlvergent flow. In this section we isolate 
the areas In which the Fourth-Order model Is 
different from the CLAS model, discuss some 
of the alternatives, and the justification of 
our choice. 

2,1 Accuracy 

Thera la a consensus among modelers that for 
finite difference models with second-order 
accuracy, horizontal resolution of about 
400 km, and about 10 vertical levels, horl- 
tontal truncation errors are the most Important 
source of errors. The truncation errors 
can ba reduced by any of the following methods: 

a. Increased horizontal resolution , 
retaining second-order differences. 

Advantages ; If the horizontal grid size 
is reduced by a factor of 2, cruncatlon 
error* are reduced by a factor of 4. A 
comparison of the effect of truncation errors 
on the computational phase speed and group 
velocity of a linear wave is presented In 
Fig. 1. Another advantage la that Increasing 
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tha resolution allow* smaller but poealbly 
important scales to be explicitly Included 
In the model. 

Disadvantages : ITte reduction of error la 

alow, and the coaputatlon tlae la lneraaaed 
by a factor of 8. 


COMPUTATIONAL PHASE I PECO efo 
10| 



I u 


i MMVntNOTW (AXI 
1 l/AX 


COMPUTATIONAL GROUP VELOCITY c'^/to 



WAVI1MOTM (AX I 
l/AX 


Fig. 1 Computational phase speed and group 
velocity fot a linear wave of 
constant phase speed c, for different 
numerical schemes as a function of 
wavelength L. 


b. Fourth order schemes of the first 
kind . We Include in this group the common 
5-polnt explicit fourth order schemes 
jKreiss and Ollger, 1972]. 

Advantages : Phase speed errors are auch 
smaller, especially for long waves (Fig. 1). 
Waves longer than 5Ax have smaller errors in 
their phase apeed with fourth order differences 
than with second order differences with twice 
the resolution. If we consider the compute* 
tlonal group velocity, the crossover occurs 
at about 7 Ax* The computational tlae la only 
Increased by a factor between 2 and 3. 
Disadvantages ; Waves In the range of 2 Ax to 
4 Ax are still grossly misrepresented- 

c. Fourth order schemes of the second 
kind. In this group we Include linear finite 
elements, cubic splines end the "compact" or 
Pade type 3-polnt 4th order difference! 
introduced by Krelss (Orstag 4 Israeli, 

1972], ail of which have the aaae fourth 
order accuracy. 


Advantagea ; These fourth-order differences 
are considerably more accurate than the 
fourth-order difference!) of the first kind 
(Fig. J). The crossover with double- resolut ion 
second-order differences in the phase speed 
error occurs at 3Ax, and in the group velocity 
error at about 4 Ax. 

Disadvantages : In the simplest case, Krelss* 

fourth order differences require the solution 
of trl-dlagonel matrices. Finite element 
schemes require the solution of at least block 
trl-dlagonal matrices. Even though there are 
efficient methods to perform these inversions, 
they are still computet tonal ly expensive. 

The extra accuracy can be* compensated with 
4th order differences of the first kind by 
Increasing the resolution, which has othor 
advantages, as we mentioned earlier. These 
schemes, ss well as in higher order finite 
difference schemes, waves In the range of 2 Ax 
to 4 Ax are still grossly misrepresented- 

d > Spectral schemes . 

Advantages : .if the basis of the spectral 
expansion are the eigenfunctions of the slave 
equation spectral sc hemes have no phase 
speed errors. Because of this they require 
less resolution than finite difference schemes. 
Disadvantages : Because of the iarge number 

of computations required for the nonlinear 
terms, spectral schemes are competitive with 
finite difference schemes only In combination 
with the use of less resolution and semi- 
Impllclt time schemes. 

Our choice : Wr chose to use fourth order 

finite differences of the first kind because 
they are computationally efficient and have 
small truncation errors except In the range 
of waves with wavelengths between 2Ax and 4Ax- 

2.2 Type of Grid 

Both staggered and unstaggered grids have been 
widely used by atmospheric modelers. 

n. Unstaggered grid : The advantage 

of this grid Is its simplicity. Higher 
order schemes are easily developed with this 
grid. Its disadvantage 1 b that all centered 
differences have to be computed over a 
distance of 2Ax. 

b. Staggered grids : Several staggered 

grid configurations are possible as reviewed 
by Arakawa [1972). The one he called scheme C, 
which Is the most commonly used, has the 
pressure defined at the center of a grid 
cell, and the velocity components u and v 
defined st their corresponding normal walls. 
This grid has the advantage that the pressure 
gradient and velocity divergence terms are 
computed over a distance of only lAx, so 
that inertia gravity vaveo are computed 
ivith double resolution. Therefore, as 
pointed out by Arakawa, geostrophlc adjust- 
ment Is best represented In this grid. On 
the other hand, advection terms are com- 
puted with no more accuracy than In the 
unstaggered grid, and In the Coriolis' 
acceleration term, it Is necessary to take 
horlcontal averagea of Che velocities. The 
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hlghtr reeolutlon of Inertia gravity wavs* 
reduces the maximum time etep for explicit 
time eohmee by a footer of two. Pull fourth 
order aeheaaa can be developed with ataggerad 
achaaaa but they are vary Involved (Kalnay- 
lUvae, 197b). To data, aodalara ualng 
ataggerad grtda have Introduced fourth order 
dlffarancaa only in the advactlon tame. 

Otir eholca i During eatended rang* fore- 
caete, aecond order arrora In the non-advectlve 
Caraa aay become important. For thla rnaaoo 
wa ehoaa to uaa an unataggarad grid that 
allova the uee of a simple, full fourth-e'rder 
achaae* 

2.3 Conservation Propertlaa an d the Dae of 

Horizontal Diffusion 

With raepect to conaervatlon propertlea. there 
are baelcally 3 type* of finite difference 
acheaea for the primitive equations: 
a) Monconserve tive schemes, the simplest of 
which is the one based on the edvectlve form 
of Che equations; b) quadrat lcally or energy 
conserving schemes; and c) enstrophy conserving 
schemes. Advectlve and quadrstically conserva- 
tive schemes can be easily developed using 
staggered or unstaggered grids (Lilly, 196S; 
Bryan, 1966). Enstrophy conserving schemes 
for the primitive equations have been developed 
on a grid C by Grammeltvedt (1969) and Arafcawa 
and Mints (1974). Sadourny (1965a, b) con- 
structed a potential enstrophy conserving 
scheme on grid C, and Arakawa (1976) has 
recently developed a potential enstrophy and 
energy conserving scheme also on grid C. 

Nonconservative schemes require a procedure 
to damp waves shorter than 4 Ax. which other- 
wise grow spuriously causing catastrophic 
nonlinear instability (Phillips, 1959). 

This has usually been do^e by means of linear 
or nonlinear horizontal diffusion, or by 
ualng dissipative numerical schemes such as 
the Lax-Wendroff , or schemes that contain 
explicit horizontal averaging. 

Quadratically conservative schemes avoid the 
unbounded growth of the solutions aasoclated 
with catastrophic nonlinear Instability. 
However, as Arakawa (1972) and Sadourny 
(1973a, b) have pointed out, in the course 
of long integrations, there Is still a spurious 
build up of energy in the shortest waves, 
which appears as an unbounded growth of the 
total enetxvphy. In the absence of hori- 
zontal diffusion, this type of alow non- 
linear instability will completely distort 
the solution. Enstrophy conserving schemes, 
on the other hsnd, impose a stronger con- 
straint on the growth of the smallest scales 
present in the model. For this reason, the 
UCLA and the GLAS models, which use enstrophy 
conserving schemes, (&> not need to include 
horieontal diffusion. 

It should be emphasized that conservation of 
enstrophy doea not necessarily Imply a more 
accurate or realistic simulation. In the 
real atmosphere, the constraint of quasl- 
^eos trophic motion Implies that very little 
of the energy generated In the barocllnlcslly 
unstsble scales can reach the smallest scales 


and ba avsntually dlaalpated (Charney, 1972). 

In a numerical modal, the finite resolution 
imposes an artificial "wall" at the short end 
of the epeotrum. Inducing an axcoaoive 
accumulation of anergy in the ahortoat wavsa. 
This problem la worst in nonconsarvstiva 
schsass, but it appeara even in alias-free, 
energy- end enetrophy-coneerving spectral 
modala. Thla Justifies ualng aoma para- 
mstarlsatlon of the unreeolved eubgrld 
eddlea to withdraw anargy from tho amalloat 
raaolvad acalea. 

For thla purpoae, Leith (1972) euggeeted the 
uee of nonlinear horizontal diffusion in 
which the addy dlffualon coefficient is pro- 
portional to the local gradient of vortlclty. 
Thla formulation la consistent with the crane- 
far of onergy to subgrid scales in two- 
dimensional turbulence. In another widely 
ueed formulation, suggested by Smagorlnsky 
(1963) and baaed on e three-dimensional 
turbulent cascade theory, the diffusion 
coefficient is proportional to the deformation 
tsnsor. These formulations are better than 
the uee of linear diffusion, bue they both 
share the following problems: a) the diffu- 

sion coefficient is computed Inaccurately 
for the shortest waves, and, more Importantly, 
b) when the diffusion coefficient is large 
enough to avoid the spurious growth of the 
smallest scales, it produces excessive 
damping of Che larger scales (Merilcea, 

1975; Williamson, 1978). Furthermore, at the 
short end of the spectrum (scales of the order 
of 100 km), neither quasi-geostrophlc nor 3- 
dimensional Isotropic turbulence theories 
are really Justified. Williamson (1978) 
generalized a higher order diffusion of the 
form , suggested by Kreiss and Oliger 

(1972:, still using the deformation type of 
diffusion coefficient. This formulation 
has the advantage that, because It is more 
scale dependent, there la less damping of 
the longer waves. 

In our model, we have taken an approach closer 
to Fourier filtering the shortest waves, 
as suggested by Phillips (1959). Our "subgrld 
parameterization’’ is based on the following 
argument: The 4th order scheme is adequate 

accurate for waves longer than 4 A* but 
grossly inaccurate for waves between 2 Ax and 
3 Ax. Since the shoreest waves cannot provide 
any ueeful information in a finite difference 
scheme, we filter them out of the system 
while their amplitude is still small. Even 
though they are filtered out, the shortest 
waves still play an Important role in the 
model: they act as a buffer or "sponge layer” 

in the spectrum domain, allowing energy eo 
trickle down from longer waves and avoiding 
tha accumulation of energy that would other- 
wlss occur at the short wave cutoff. 

The elimination of the short waves is performed 
In the model with the periodic use of a l 6th 
order Shapiro (1970) filter, first Introduced 
in GCMa by Francis (1975). Figure 2 indi- 
cates the response of Shapiro filters of 
order 4, 8 and 16. It can be observed that 
in the case of a 16th order filter waves 
longer than 4 are scarcely affected at all 
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•ven after 128 applications, which In our 
■odel correspond to a 10 day integration. 
Waves shorter than 4 are virtually ellal- 
nated. Lower order filters Introduce too 
■uch decay at long scales. 





Fig. 2 Response of Shapiro filters after 1 
and after 128 applications. 


In order to explore the effect of filtering 
short waves we performed a series of 20 day 
forecasts with a simple shallow watar equa- 
tion model [Kalnay-Rlvas, 19761, using diffe- 
rent combinations of schemes, smoothing 
operators and frequency of application as 
Indicated on Table 1. 


Alt 

4th order. 

non-conservative, no smoothing. 

Alt 

•9 

• 

•I 

.linear dlf fusion •O Ws. 

A3: 

if 

B 

M 

,16th order filter/4 hrs. 

A4: 

M 

• 

•• 

, " /time step. 

Bit 

M 

» ' 

quadrat 

« conaerv., no smoothing. 

B2; 

M 

t 

If 

,16th order fllter/4 hrs. 

B3: 

if 

B 

If 

, 8th " 

B4; 

*# 

t 

It 

, 4th * 

C: 

2nd ordmr. 

91 

, no smoothing 


TABLE 1: Characteristic* of the different 
runs In Fig. 3. 
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Fig. 3 Variations of total enstrophy during 
numerical integrations as indicated 
In Table 1. The scale corresponding 
to experiments A1 and A2 Is multi- 
plied by 10“ . 

We found that all stable runs conserved total 
energy with a high degree of accuracy. Fig. 3 
shows the variation lr time of the total 
potential enstrophy , which Is conserved 
exactly In the continuous equations. The 
results Indicate that when the 16th order 
filter is applied every time step (10 minutes) 
even a formally nonconsevvatlve scheme con- 
serves both total energy and total potential 
enstrophy during a long Integration (Run A4). 
The quadratlcally conservative scheme controls 
better the amount of energy going Into the 
smallest scales, so that in Run B2 It was 
enough to apply the filter every 4 hours to 
conserve potential enstrophy within 0.052, 
even though such conservation is hot formally 
guaranteed in the scheme. 
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It any bo questioned whether the application 
of a high order filter eliminates swell 
scale features like frontal zones or the 
effects of orographic or cumulus convection 
forcing on snail scales. Pig. 4 Indicates 
that this Is true for lover order filters* 
However a 16th order filter ellnlnatee only 
those components which are not resolved any- 
way, and still allows for the forest Ion of 
sharp gradient cones end strong local sax lea. 

10 PASSES (ONE DAY) OP SHAPIRO FILTER 


Our choice t We choee to use e quadrat* 
ically (energy) conserving scheme because It 
la only slightly less efficient then an advec- 
tlva scheme end It requires the application 
of a hlgh-order filter only every few houra 
to avoid the effect of slow nonlinear Insta- 
bility associated with a spurious growth 
of snatrophy* Wa apply the filter to the 
gga level preaeure and potential tempera- 
ture fielde In order to compensate for the 
effect of topography. 




3. EXPERIMENTAL FORECASTS: PRELIMINARY 

RESULT S 

We plan to perfors an extensive series of 
forecasts to study the effect of using 
different schemes and varying resolution on 
tha quality of aetual wether forecasts. In 
this section we present suae preliminary 
results. 

We tested the Fourth Order model by making 3- 
day forecasts from several initial conditions. 
In every case the model performed much better 
then the CLAS model with the same 4* by 3* 
resolution. After 3 days the forecasts were 
comparable to the 2.3* by 3* version of the 
CLAS model in the sea level pressure maps, 
and had slightly less phase errors in the 
300 ab maps. 

Fig. Sa shows the verification sea level 
pressure map corresponding to a 3-day forecast 
with February 19, 1976, 0Z (a case that has 
been studied In detail by Atlas et al. . 1979). 
Fig. Sb Is the 3-day forecast with the Fourth 
Order 4* by 5* model verifying on the same 
date. The excessive gradients, especially at 
high latitudes are due to a faulty computation 
of the ground temperature by the radiation 
routine used in that run. Figs. 3c and 5d 
are the forecasts generated by the CLAS models 
with 4* by S* and 2.5* by 3* resolution res- 
pectively. Fig. 6 displays the corresponding 
300 ab maps. 

In Fig. 7a, we present the verification sea 
level pressure map corresponding to a 3-day 
forecast with February 1, 1976, 0Z initial 
conditions. Hie results of six different 3- 
day forecasts are shown in Figs. 7b to 7g. 

The forecast in Fig. 7b was computed with the 
new full Fourth"Order model, using the same 
resolution and parameterization of physical 
processes (except for the long-wave radiation 
routine) as in the 4* by 5* GLAS model 
(Fig. 7c). In Fig. 7d, the model was the 
same as in 7b, but full Second-Order accuracy 
wee used. . If we compare these three made 
with the same resolution, we see that the 
fourth order model did considerably better, 
especially In forecasting the development and 
aotlon of the cyclone southwest of Greenland. 
The two second-order foreceats are close to 
each other. 


Fig. 4 


<3> 


Effect of 10 applications of a 16th 
order Shapiro filter on a) a step 
function b) a spike. 
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The idee of filtering these fielde is due 
to Dr. A. Bayliss. 



Fig. 7 shows the foraeaat mad* with tha 2«S* 
by 3* GLAS model atardng from the NASA 
initial conditions with aaalallation of 
satellite data (Atlaa ft al. , 1979). All 
other forecaata were made from NMC'a Global 
Analysis Initial conditions. The cyclo- 
genesis forecast was poorer than with tha 
fourth order model, but over continental 
North America the forecaat was better. This 
stay be due to the higher resolution or, 
possibly, Improved Initial conditions. 

Fig. 7f presents the forecast made with NMC’a 
6-layer, 380 km resolution, which has resolution 
comparable to our 4* by 5* grid. For theae 
initial conditions, NMC’a forecaat waa better 
than the GLAS second order forecast, although 
it snared some of Its errors (such as a 
spurious antlcyclogenesla over the Great 
Lakes). 

Fig. 7g shows a forecast made with a slightly 
different version of the fourth order model. 

The differences were as follows; In the 7g 
forecast we used a Matsuno time step, a 
surface drag that was .75 of that of the 
standard model, the Shapiro was applied every 
two hours, and a slightly different scheme 
for the vertical advectlon of moisture wee 
used. In the 7b forecast we used a combined 
Mateuno-leap-f rog scheme, the surface drag 
was the same as In the standard model and the 
filter was applied every hour. The two 
forecasts are quite similar, but from other 
experiments it seems that the positions of 
the oceanic lows west of Spain and south of 
Alaska were somewhat affected by the change 
In frequency of the Shapiro filter. 


4. FINAL REMARKS 

Although these are preliminary results, the 
new Fourth-Order model aeems to have very 
good forecasting skill. Similar excellent 
results with fourth order* schemes were reported 
by Camparia (1978, 1979), and by Williamson 
(1978). Campana found that most of the Improve- 
ment over the second order model was obtained 
Just by Introducing fourth order accuracy on 
the horlcontal averages performed on the 
Shuman-Hove rmftle model. Campana also obtained 
that for 2-day forecasts, fourth order 
differences were important In the advectlon 
terms but not In the pressure and continuity 
terms. We made shallow water equation 
experiments that showed some deterioration 
of the solution after about 8 days when only 
the advectlon terms wera written with fourth 
order accuracy. Since barocllnlc models are 
much more unstable than the SWE, Campana' s 
results may not hold during extended forecasts. 
We are performing experiments to study this 
possibility In our model. We are also 
repeating some of the forecasts using a 
vector Invariant form of the momentum equa- 
tions. The NMC 6-layer PE model used In 
Fig. 7f Is based on such a formulation, and 
we want to determine If It affects the 
solution. 


We are performing Integrations with a 2.5* by 
3* version of the Fourth-Order model. From 
linear theory, this should make horlcontal 
truncation errors negligible for most synoptic 
scales for periods of a week or longer* 

Still, narrow atmospheric features like those 
due to sharp orographic forcing on the ITCZ 
will remain poorly resolved. 

We have found the fourth order model to ba 
extremely sensitive to the parsmeterlcstlon 
of physical processes. For example. It beeame 
unstable through excessive cooling near the 
surface, apparently due to a flaw In the 
radiation routine. None of the second order 
models was sensitive to this problem. 

In this paper we have discussed the use of 
formally nonconservatlve schemes coupled with 
periodic filtering of the shortest waves as 
an alternative to the use of conservative 
schemes. It Is clear from the experimental 
results that formal enstrophy conservation 
has had no beneficial Impact on the quality 
of the forecast (Figs. 7c snd 7d). Conservation 
of potential enstrophy (Arekawa, 1978] might 
have a more positive effect. 

We think that the use of nonconservatlve 
schemes with high order periodic filtering le 
also Justified In climate simulations as long 
as the rates of energy or enstrophy loss due 
to filtering remain much smaller than the 
observed rates of generation and dissipation. 
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Fig. 5 Sea level pressure maps corres- 
ponding to 3-day forecasts with 
February 19, 1976, 02 Initial 
conditions* 

Fig. 6 Same ss Fig. 5 but 500 mb 
geopotsntlsl heights. 

Fig. 7 Same ss Fig. 5 but for February 1, 
1976, OZ Initial conditions. 
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